OCS ePS results processing + AF-BLMs (Jan 2022 update)

Quick attempt at AF-BLMs for OCS, based on guessed/reasonable ADMs - action starts below after setup stage). Status: basically working, may be some remaining bugs with ADMs, cos^2 and renormalisation, but trends with alignment should be correct.

For methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html

Additional notes below

Setup

In [1]:
!hostname
jake
In [2]:
!conda env list
# conda environments:
#
base                     /home/paul/anaconda3
automation               /home/paul/anaconda3/envs/automation
baseBK-250821            /home/paul/anaconda3/envs/baseBK-250821
baseBK-280321            /home/paul/anaconda3/envs/baseBK-280321
ePSproc-v1.2             /home/paul/anaconda3/envs/ePSproc-v1.2
ePSproc-v1.2-dev         /home/paul/anaconda3/envs/ePSproc-v1.2-dev
epsTestConda             /home/paul/anaconda3/envs/epsTestConda
epsTestPip               /home/paul/anaconda3/envs/epsTestPip
epsTestPypi              /home/paul/anaconda3/envs/epsTestPypi
epsTestSetup             /home/paul/anaconda3/envs/epsTestSetup
epsdev                   /home/paul/anaconda3/envs/epsdev
epsdev-040821            /home/paul/anaconda3/envs/epsdev-040821
epsdev-040821-YMLnb      /home/paul/anaconda3/envs/epsdev-040821-YMLnb
epsdev-xr13              /home/paul/anaconda3/envs/epsdev-xr13
epsdev-xr15           *  /home/paul/anaconda3/envs/epsdev-xr15
epsdev-xr17              /home/paul/anaconda3/envs/epsdev-xr17
epsdevTestYML            /home/paul/anaconda3/envs/epsdevTestYML
epsman                   /home/paul/anaconda3/envs/epsman
fibre-sim                /home/paul/anaconda3/envs/fibre-sim
frog                     /home/paul/anaconda3/envs/frog
matlab                   /home/paul/anaconda3/envs/matlab
tmo-dev                  /home/paul/anaconda3/envs/tmo-dev

In [3]:
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr

import matplotlib.pyplot as plt

from datetime import datetime as dt
timeString = dt.now()

import epsproc as ep

# Plotters
import hvplot.xarray
from epsproc.plot import hvPlotters

# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob

plotWidth = 700
hvPlotters.setPlotters(width = plotWidth, snsStyle='whitegrid')
hvPlotters.opts.defaults(hvPlotters.opts.Curve(height=plotWidth))   # Fix plot height! Currently not set by default.
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* pyevtk not found, VTK export not available. 
In [4]:
# For class, above settings don't take, not sure why, something to do with namespaces/calling sequence?
# Overriding snsStyle does work however... although NOT CONSISTENTLY????
# AH, looks like ordering matters - set_style LAST (.set seems to override)
import seaborn as sns

sns.set(rc={'figure.figsize':(10,6)})  # Set figure size in inches
sns.set_context("paper")
sns.set_style("whitegrid")  # Set plot style
sns.set_palette("Paired")   # Set colour mapping

# Try direct fig type setting for PDF output figs
from IPython.display import set_matplotlib_formats
# set_matplotlib_formats('png', 'pdf')
set_matplotlib_formats('svg', 'pdf')
/home/paul/anaconda3/envs/epsdev-xr15/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)
In [5]:
# xr.set_options(display_style='html')

Load data

In [6]:
import warnings
# warnings.filterwarnings('once')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
warnings.filterwarnings('ignore')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
In [7]:
# # Scan for subdirs, based on existing routine in getFiles()

fileBase = Path('/home/paul/ePS/OCS/OCS_survey')  # Data dir on Stimpy
In [8]:
# TODO: fix orb label here, currently relies on (different) fixed format

data = ePSmultiJob(fileBase, verbose = 0)

data.scanFiles()
data.jobsSummary()
*** Warning: Missing records, expected 64, found 48.
*** Warning: Found 16 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 64, found 48.
*** Warning: Found 16 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 64, found 48.
*** Warning: Found 16 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 64, found 48.
*** Warning: Found 16 blank sets of matrix elements, symmetries ['A2']
Found 4 directories, with 8 files.

*** Job orb11 details
Key: orb11
Dir /home/paul/ePS/OCS/OCS_survey/orb11, 2 file(s).
{   'batch': 'ePS OCS, batch OCS_survey, orbital orb11',
    'event': ' orb 11 ionization, basic survey run.',
    'orbE': -11.35803261907539,
    'orbLabel': '# OCS, orb 11 ionization, basic survey run.'}

*** Job orb9 details
Key: orb9
Dir /home/paul/ePS/OCS/OCS_survey/orb9, 2 file(s).
{   'batch': 'ePS OCS, batch OCS_survey, orbital orb9',
    'event': 'orb 9 (E1/P) ionization, basic survey run.',
    'orbE': -18.34863774566971,
    'orbLabel': 'E1/P'}

*** Job orb10 details
Key: orb10
Dir /home/paul/ePS/OCS/OCS_survey/orb10, 2 file(s).
{   'batch': 'ePS OCS, batch OCS_survey, orbital orb10',
    'event': 'orb 10 (A1/S) ionization, basic survey run.',
    'orbE': -17.32548962282056,
    'orbLabel': 'A1/S'}

*** Job orb8 details
Key: orb8
Dir /home/paul/ePS/OCS/OCS_survey/orb8, 2 file(s).
{   'batch': 'ePS OCS, batch OCS_survey, orbital orb8',
    'event': 'orb 8 (A1/S) ionization, basic survey run.',
    'orbE': -21.51332196607811,
    'orbLabel': 'A1/S'}
In [9]:
# Fix labels for plots - currently have some mis-named files!

# Set state labels
v = ['C','B','A','X']
sDict = {n:v[k] for k,n in enumerate(range(9,13))}

for key in data.data.keys():
    data.data[key]['jobNotes']['orbKey'] = key  # Existing orb key, from filename
    data.data[key]['jobNotes']['orbGroup'] = int(key.strip('orb')) + 1   # Corrected orb group
#     data.data[key]['jobNotes']['orbLabel'] = f"HOMO - {12 - int(key.strip('orb')) - 1}"
    data.data[key]['jobNotes']['orbLabel'] = sDict[data.data[key]['jobNotes']['orbGroup']]
    
#     print(data.data[key]['jobNotes']['orbLabel'])

System properties

Note orbital numbering in table below:

  • Orb is Gamess file output orbital numbering, but energy but not grouped by degeneracy.
  • OrbGrp is grouped numbering by degeneracy, also used by ePolyScat, and will be used for labels etc. below.

BUT - file names are 0-indexed (oops), so give OrbGrp-1 here. Sorry.

In [10]:
data.molSummary()
*** Molecular structure
2022-01-07T13:46:45.562382 image/svg+xml Matplotlib v3.3.4, https://matplotlib.org/
*** Molecular orbital list (from ePS output file)
EH = Energy (Hartrees), E = Energy (eV), NOrbGrp, OrbGrp, GrpDegen = degeneracies and corresponding orbital numbering by group in ePS, NormInt = single centre expansion convergence (should be ~1.0).
props Sym EH Occ E NOrbGrp OrbGrp GrpDegen NormInt
orb
1 S -91.9901 2.0 -2503.178142 1.0 1.0 1.0 0.660569
2 S -20.6748 2.0 -562.589968 1.0 2.0 1.0 0.963189
3 S -11.4451 2.0 -311.437037 1.0 3.0 1.0 0.999998
4 S -8.9936 2.0 -244.728323 1.0 4.0 1.0 0.951187
5 S -6.6764 2.0 -181.674099 1.0 5.0 1.0 0.993503
6 P -6.6721 2.0 -181.557090 1.0 6.0 2.0 0.978550
7 P -6.6721 2.0 -181.557090 2.0 6.0 2.0 0.978550
8 S -1.5373 2.0 -41.832064 1.0 7.0 1.0 0.998241
9 S -1.0876 2.0 -29.595104 1.0 8.0 1.0 0.997329
10 S -0.7906 2.0 -21.513322 1.0 9.0 1.0 0.999269
11 P -0.6743 2.0 -18.348638 1.0 10.0 2.0 0.999887
12 P -0.6743 2.0 -18.348638 2.0 10.0 2.0 0.999887
13 S -0.6367 2.0 -17.325490 1.0 11.0 1.0 0.998566
14 P -0.4174 2.0 -11.358033 1.0 12.0 2.0 0.998839
15 P -0.4174 2.0 -11.358033 2.0 12.0 2.0 0.998839
*** Warning: some orbital convergences outside single-center expansion convergence tolerance (0.01):
[[1.         0.66056934]
 [2.         0.96318872]
 [4.         0.95118653]
 [6.         0.9785498 ]
 [7.         0.9785498 ]]
In [11]:
# Fix orb names! NOW FIXED ABOVE
/home/paul/anaconda3/envs/epsdev-xr15/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)

AF-$\beta_{LM}$ (UPDATES 13/12/21)

13/12/21 PH

Revisiting AF cases with (hopefully) better match to experimental alignment.

  • Add in cos^2(theta) metric.
  • Ballpark ADMs assuming Gaussian case.
  • Calculate & plot AF-BLMs for these cases.

Alignment metrics

Here we'll start with the molecular alignment defined in terms of axis distribution moments (ADMs), given by parameters $A^K_{Q,S}(t)$:

$ P(\Omega,t) = \sum_{K,Q,S} A^K_{Q,S}(t)D^K_{Q,S}(\Omega)$

Where $P(\Omega,t)$ is the full axis distribution probability, expanded for Euler angles $\Omega$.

This reduces to the 2D case if $S=0$:

$ P(\theta,\phi,t) = \sum_{K,Q} A^K_{Q,0}(t)D^K_{Q,0}(\Omega) = \sum_{K,Q} A^K_{Q}(t)Y_{K,Q}(\Omega)$

In the current code (v1.3.0, 16/08/21), the ADMs can be set directly, and expanded trivially for the 2D case (3D case to do). (Note that the 3D case is supported for AF-BLM calculations, but not yet for direct expansion & visualisation.)

Expectation values of $\langle cos^n(\theta, t)\rangle$ can be calculated directly from the $P(\theta,t)$ distributions:

$\langle cos^n(\theta, t)\rangle = \int_{\theta} cos^n(\theta)P(\theta,t) d\theta$

See, e.g.:

[1]
P. Hockett, “General phenomenology of ionization from aligned molecular ensembles,” New Journal of Physics, vol. 17, no. 2, p. 023069, Feb. 2015, doi: 10.1088/1367-2630/17/2/023069.
[2] Reid, K.L. (2018) ‘Accessing the molecular frame through strong-field alignment of distributions of gas phase molecules’, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 376(2115), p. 20170158. doi:10.1098/rsta.2017.0158.

Set test ADMs

Roughly eyeballed for a Gaussian distribution with $0<\sigma^{2}<0.2$ as per Figure 1 in Reid 2018 [2].

NOTE: not totally sure on the normalization here. Assuming othonorm sph, may be an additional sqrt(4pi) required for Y00? https://shtools.github.io/SHTOOLS/complex-spherical-harmonics.html#orthonormalized

UPDATE: setting for 4$pi$- normalised & converting to orthonorm (for numerical functions) looks good for comparison with ref. [2]. Note also additional factor of sqrt(2) in metrics at the moment - due to 0:$pi$ default domain and/or renorm(?)... should also confirm there is no additional ^2 in distribution generation.

In [12]:
# Set ADMs for increasing alignment...
tPoints = 10
# inputADMs = np.asarray([[0,0,0, *np.ones(10)], [2,0,0, *np.linspace(0,1,tPoints)], [4,0,0, *np.linspace(0,0.5,tPoints)]])
# inputADMs = [[0,0,0, *np.ones(10) * np.sqrt(4*np.pi)], 
#              [2,0,0, *np.linspace(1,2,tPoints)], 
#              [4,0,0, *np.linspace(0.3,2.8,tPoints)],
#              [6,0,0, *np.linspace(0,3,tPoints)],
#              [8,0,0, *np.linspace(0,3,tPoints)]]

# Ballpark ADMs from Fig. 1 in [2]
inputADMs = [[0,0,0, *np.ones(10) * np.sqrt(4*np.pi)], 
             [2,0,0, *np.linspace(0.3,2.2,tPoints)], 
             [4,0,0, *np.linspace(0,1.5,tPoints)],
             [6,0,0, *np.linspace(0,1.5,tPoints)],
             [8,0,0, *np.linspace(0,1.5,tPoints)]]


warnings.filterwarnings('ignore')   # Skip repeated numpy deprecation warnings in current build (xr15 env)

ADMs = ep.setADMs(ADMs = inputADMs)  # TODO WRAP TO CLASS (IF NOT ALREADY!)

# Set to separate dataset for testing
ADMs.attrs['jobLabel'] = 'Alignment data testing'
data.data['ADM'] = {'ADM': ADMs}

ep.lmPlot(ADMs, xDim='t');
Plotting data (No filename), pType=a, thres=0.01, with Seaborn
No handles with labels found to put in legend.
2022-01-07T13:46:46.024803 image/svg+xml Matplotlib v3.3.4, https://matplotlib.org/

Compute $P(\theta)$ distribtuions + metrics

In [13]:
# Axis distributions P(theta) vs. t (summed over phi)
# Pt = data.padPlot(keys = 'ADM', dataType='ADM', Etype = 't', pStyle='grid', reducePhi='sum', returnFlag = True)  #, facetDims=['t', 'Eke'])
Pt = data.padPlot(keys = 'ADM', dataType='ADM', Etype = 't', pType='r', pStyle='grid', reducePhi='sum', returnFlag = True)  #, facetDims=['t', 'Eke'])  # Real plot to check values OK.
Using default sph betas.
Summing over dims: set()
Grid plot: Alignment data testing, dataType: ADM, plotType: r
2022-01-07T13:46:46.469943 image/svg+xml Matplotlib v3.3.4, https://matplotlib.org/
In [14]:
# The returned data array contains the plotted values
Pt
Out[14]:
Show/Hide data repr Show/Hide attributes
xarray.DataArray
  • t: 10
  • Theta: 50
  • 59.46 59.4 59.23 58.94 58.55 58.05 ... 199.4 255.3 303.0 335.1 346.4
    array([[ 59.46174696,  59.40348635,  59.22966116,  58.9431256 ,
             58.54858456,  58.05251641,  57.46306657,  56.78991379,
             56.04411122,  55.23790491,  54.38453276,  53.49800711,
             52.59288468,  51.68402756,  50.78635916,  49.91461917,
             49.08312154,  48.30551944,  47.59458106,  46.96198   ,
             46.41810356,  45.97188216,  45.63064274,  45.39998846,
             45.28370665,  45.28370665,  45.39998846,  45.63064274,
             45.97188216,  46.41810356,  46.96198   ,  47.59458106,
             48.30551944,  49.08312154,  49.91461917,  50.78635916,
             51.68402756,  52.59288468,  53.49800711,  54.38453276,
             55.23790491,  56.04411122,  56.78991379,  57.46306657,
             58.05251641,  58.54858456,  58.9431256 ,  59.22966116,
             59.40348635,  59.46174696],
           [ 91.34083196,  90.03125136,  86.309706  ,  80.75509747,
             74.19972596,  67.56228639,  61.669788  ,  57.10773617,
             54.13034905,  52.64770197,  52.28823463,  52.51757093,
             52.78231608,  52.64326306,  51.86704996,  50.45732705,
             48.62266126,  46.69451041,  45.02055356,  43.86355124,
             43.33271181,  43.36435941,  43.75443271,  44.23090009,
             44.54343007,  44.54343007,  44.23090009,  43.75443271,
             43.36435941,  43.33271181,  43.86355124,  45.02055356,
             46.69451041,  48.62266126,  50.45732705,  51.86704996,
             52.64326306,  52.78231608,  52.51757093,  52.28823463,
             52.64770197,  54.13034905,  57.10773617,  61.669788  ,
             67.56228639,  74.19972596,  80.75509747,  86.309706  ,
             90.03125136,  91.34083196],
           [123.21991695, 120.65901637, 113.38975085, 102.56706934,
             89.85086735,  77.07205637,  65.87650943,  57.42555855,
             52.21658688,  50.05749903,  50.19193651,  51.53713476,
             52.97174749,  53.60249856,  52.94774076,  51.00003493,
             48.16220097,  45.08350138,  42.44652607,  40.76512248,
             40.24732007,  40.75683666,  41.87822267,  43.06181172,
             43.80315349,  43.80315349,  43.06181172,  41.87822267,
             40.75683666,  40.24732007,  40.76512248,  42.44652607,
             45.08350138,  48.16220097,  51.00003493,  52.94774076,
             53.60249856,  52.97174749,  51.53713476,  50.19193651,
             50.05749903,  52.21658688,  57.42555855,  65.87650943,
             77.07205637,  89.85086735, 102.56706934, 113.38975085,
            120.65901637, 123.21991695],
           [155.09900195, 151.28678137, 140.46979569, 124.37904122,
            105.50200875,  86.58182635,  70.08323086,  57.74338094,
             50.30282471,  47.46729609,  48.09563839,  50.55669858,
             53.16117889,  54.56173405,  54.02843155,  51.54274281,
             47.70174068,  43.47249235,  39.87249857,  37.66669371,
             37.16192832,  38.14931392,  40.00201263,  41.89272336,
             43.06287692,  43.06287692,  41.89272336,  40.00201263,
             38.14931392,  37.16192832,  37.66669371,  39.87249857,
             43.47249235,  47.70174068,  51.54274281,  54.02843155,
             54.56173405,  53.16117889,  50.55669858,  48.09563839,
             47.46729609,  50.30282471,  57.74338094,  70.08323086,
             86.58182635, 105.50200875, 124.37904122, 140.46979569,
            151.28678137, 155.09900195],
           [186.97808695, 181.91454638, 167.54984054, 146.19101309,
            121.15315014,  96.09159632,  74.28995229,  58.06120332,
             48.38906254,  44.87709314,  45.99934026,  49.57626241,
             53.3506103 ,  55.52096955,  55.10912235,  52.08545069,
             47.24128039,  41.86148331,  37.29847107,  34.56826495,
             34.07653658,  35.54179117,  38.12580259,  40.72363499,
             42.32260034,  42.32260034,  40.72363499,  38.12580259,
             35.54179117,  34.07653658,  34.56826495,  37.29847107,
             41.86148331,  47.24128039,  52.08545069,  55.10912235,
             55.52096955,  53.3506103 ,  49.57626241,  45.99934026,
             44.87709314,  48.38906254,  58.06120332,  74.28995229,
             96.09159632, 121.15315014, 146.19101309, 167.54984054,
            181.91454638, 186.97808695],
           [218.85717195, 212.54231139, 194.62988538, 168.00298497,
            136.80429154, 105.6013663 ,  78.49667372,  58.3790257 ,
             46.47530037,  42.2868902 ,  43.90304214,  48.59582624,
             53.5400417 ,  56.48020505,  56.18981315,  52.62815857,
             46.7808201 ,  40.25047428,  34.72444358,  31.46983618,
             30.99114483,  32.93426843,  36.24959256,  39.55454662,
             41.58232376,  41.58232376,  39.55454662,  36.24959256,
             32.93426843,  30.99114483,  31.46983618,  34.72444358,
             40.25047428,  46.7808201 ,  52.62815857,  56.18981315,
             56.48020505,  53.5400417 ,  48.59582624,  43.90304214,
             42.2868902 ,  46.47530037,  58.3790257 ,  78.49667372,
            105.6013663 , 136.80429154, 168.00298497, 194.62988538,
            212.54231139, 218.85717195],
           [250.73625695, 243.1700764 , 221.70993023, 189.81495684,
            152.45543293, 115.11113628,  82.70339515,  58.69684808,
             44.5615382 ,  39.69668726,  41.80674401,  47.61539006,
             53.72947311,  57.43944055,  57.27050394,  53.17086645,
             46.32035981,  38.63946525,  32.15041608,  28.37140742,
             27.90575309,  30.32674568,  34.37338252,  38.38545825,
             40.84204718,  40.84204718,  38.38545825,  34.37338252,
             30.32674568,  27.90575309,  28.37140742,  32.15041608,
             38.63946525,  46.32035981,  53.17086645,  57.27050394,
             57.43944055,  53.72947311,  47.61539006,  41.80674401,
             39.69668726,  44.5615382 ,  58.69684808,  82.70339515,
            115.11113628, 152.45543293, 189.81495684, 221.70993023,
            243.1700764 , 250.73625695],
           [282.61534195, 273.79784141, 248.78997507, 211.62692872,
            168.10657432, 124.62090626,  86.91011658,  59.01467046,
             42.64777603,  37.10648432,  39.71044589,  46.63495389,
             53.91890451,  58.39867605,  58.35119474,  53.71357432,
             45.85989952,  37.02845622,  29.57638858,  25.27297866,
             24.82036134,  27.71922294,  32.49717248,  37.21636988,
             40.1017706 ,  40.1017706 ,  37.21636988,  32.49717248,
             27.71922294,  24.82036134,  25.27297866,  29.57638858,
             37.02845622,  45.85989952,  53.71357432,  58.35119474,
             58.39867605,  53.91890451,  46.63495389,  39.71044589,
             37.10648432,  42.64777603,  59.01467046,  86.91011658,
            124.62090626, 168.10657432, 211.62692872, 248.78997507,
            273.79784141, 282.61534195],
           [314.49442695, 304.42560642, 275.87001992, 233.43890059,
            183.75771572, 134.13067623,  91.11683801,  59.33249284,
             40.73401386,  34.51628138,  37.61414776,  45.65451772,
             54.10833592,  59.35791155,  59.43188553,  54.2562822 ,
             45.39943924,  35.41744719,  27.00236108,  22.17454989,
             21.7349696 ,  25.11170019,  30.62096244,  36.04728151,
             39.36149402,  39.36149402,  36.04728151,  30.62096244,
             25.11170019,  21.7349696 ,  22.17454989,  27.00236108,
             35.41744719,  45.39943924,  54.2562822 ,  59.43188553,
             59.35791155,  54.10833592,  45.65451772,  37.61414776,
             34.51628138,  40.73401386,  59.33249284,  91.11683801,
            134.13067623, 183.75771572, 233.43890059, 275.87001992,
            304.42560642, 314.49442695],
           [346.37351194, 335.05337143, 302.95006476, 255.25087246,
            199.40885711, 143.64044621,  95.32355944,  59.65031522,
             38.82025169,  31.92607843,  35.51784964,  44.67408154,
             54.29776732,  60.31714704,  60.51257633,  54.79899008,
             44.93897895,  33.80643816,  24.42833359,  19.07612113,
             18.64957786,  22.50417744,  28.74475241,  34.87819314,
             38.62121744,  38.62121744,  34.87819314,  28.74475241,
             22.50417744,  18.64957786,  19.07612113,  24.42833359,
             33.80643816,  44.93897895,  54.79899008,  60.51257633,
             60.31714704,  54.29776732,  44.67408154,  35.51784964,
             31.92607843,  38.82025169,  59.65031522,  95.32355944,
            143.64044621, 199.40885711, 255.25087246, 302.95006476,
            335.05337143, 346.37351194]])
    • t
      (t)
      int64
      0 1 2 3 4 5 6 7 8 9
      units :
      Index
      array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    • Theta
      (Theta)
      float64
      0.0 0.06411 0.1282 ... 3.077 3.142
      array([0.      , 0.064114, 0.128228, 0.192342, 0.256457, 0.320571, 0.384685,
             0.448799, 0.512913, 0.577027, 0.641141, 0.705255, 0.76937 , 0.833484,
             0.897598, 0.961712, 1.025826, 1.08994 , 1.154054, 1.218169, 1.282283,
             1.346397, 1.410511, 1.474625, 1.538739, 1.602853, 1.666968, 1.731082,
             1.795196, 1.85931 , 1.923424, 1.987538, 2.051652, 2.115766, 2.179881,
             2.243995, 2.308109, 2.372223, 2.436337, 2.500451, 2.564565, 2.62868 ,
             2.692794, 2.756908, 2.821022, 2.885136, 2.94925 , 3.013364, 3.077479,
             3.141593])
  • dataType :
    ADM
    long_name :
    Axis distribution moments
    units :
    arb
    jobLabel :
    Alignment data testing
    pType :
    r
    pTypeDetails :
    {'Type': 'real', 'Formula': 'np.real(dataPlot)', 'Exec': <function plotTypeSelector.<locals>.<lambda> at 0x7f6d3c882b90>}
In [15]:
# Check no negative values...
print(Pt.min() > 0)
<xarray.DataArray ()>
array(True)
In [16]:
def computeCNs(Pt, N = None, Nrange = [2,10], Nstep = 2):
    """
    Function to compute <cos^N> expectation values from P(theta) distributions.
    
    NOTE: currently assumes 1D case. Should generalise!
    
    """
    
    if N is None:
        N = np.arange(Nrange[0], Nrange[1], Nstep)
    
    # Generate cos^N basis
    #  cn = xr.ones_like(Pt) * (np.cos(Pt.Theta)**N)  # Use multiplication to keep t coord - this fails for array-like N with broadcast errors
    cn = xr.ones_like(Pt) * np.cos(Pt.Theta).expand_dims({'N':N}).T.pipe(np.power,N)  # OK with transpose
    
    Jt = np.sin(Pt.Theta) # Jacobian
    
    norm = (Pt * Jt).sum('Theta')
    cnMetric = (Pt * cn * Jt).sum('Theta')/norm
    
#     cnMetric = Pt.dot(cn, dims='Theta') # Multiply and sum over Theta - dot version (no Jacobian or renorm)
#     cnMetric = (Pt.dot(cn * Jt, dims='Theta')/Pt.dot(Jt, dims='Theta'))  # Dot version with renorm
    
    return cnMetric
In [17]:
# Compute and plot for a range of N
PtNorm = Pt  #/Pt.sum('Theta')  # Norm integral to 1 - now set in function
cnMetric = computeCNs(PtNorm)

# cnMetric.plot()  # Basic plot
sf = np.sqrt(2)  # Seem to need additional sqrt(2) normalisation? Probably SPH definition and/or integration domain and/or distribution generation params.
# sf = 1
ep.lmPlot(cnMetric*sf, xDim='t');  # With lmPlot()
Set dataType (No dataType)
Plotting data (No filename), pType=a, thres=0.01, with Seaborn
No handles with labels found to put in legend.
2022-01-07T13:46:47.198316 image/svg+xml Matplotlib v3.3.4, https://matplotlib.org/
In [18]:
# Seem to need additional sqrt(2) normalisation? Probably SPH definition?
# (cnMetric*np.sqrt(2)).sel(N=2).plot()
# (cnMetric*np.sqrt(2)).plot.line(x='t')

# Testing exact values
# sf = np.sqrt(2)
# sf = 2
# (cnMetric*sf).plot.line(x='t')
# (ADMs).plot.line(x='t')

(cnMetric*sf).fillna(0).hvplot.line(x='t').overlay('N')
Out[18]:

AFBLMs for aligned case

New plotting code to update current package versions!

In [19]:
# Compute AFBLMs for aligned case
keys = data._keysCheck(None)
keys = list(set(keys) - set(['ADM']))

data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'})
In [20]:
# Quick stack to dict & XR dataset
dataType = 'AFBLM'
thres = 1e-2
selDims = None
Etype = 'Eke'

pDict = []

for key in keys:
    subset = ep.matEleSelector(data.data[key][dataType], thres=thres, inds = selDims, dims = Etype, sq = True)
    subset.name = 'BLM'

    subset = ep.plotTypeSelector(subset, pType = 'r')

#     pDict.append(subset.expand_dims({'Orb':[key]}))
#     pDict.append(subset.expand_dims({'Orb':[data.data[key][dataType].attrs['jobLabel']]}))
    pDict.append(subset.expand_dims({'Orb':[data.data[key]['jobNotes']['orbLabel']]}))

    
#     hvDS = hv.Dataset(subset.unstack('BLM'))
    
#     display(hvDS.to(hv.Curve, kdims=Etype).overlay(['l','m']))

#     pDict[key] = hvDS.to(hv.Curve, kdims=Etype).overlay(['l','m'])
    
# hv.HoloMap(pDict, kdims = ['Orb'])   # Individual plots OK, but won't stack?

xrAF = xr.concat(pDict, 'Orb')
In [21]:
# hvPlotters.opts.defaults(hvPlotters.opts.Curve(height=700))
hvDS = hvPlotters.hv.Dataset(xrAF.unstack('BLM'))
hvDS.select(l=[2,4,6]).to(hvPlotters.hv.Curve, kdims=Etype).overlay(['l','m'])
Out[21]:
In [22]:
# With additional overlays
hvDS.select(l=[2,4,6], m=0).to(hvPlotters.hv.Curve, kdims=Etype).overlay(['t']).layout('Orb').cols(1)
Out[22]:
In [23]:
# With additional overlays
# This seems to break legend?
hvDS.select(l=[2,4,6]).to(hvPlotters.hv.Curve, kdims=Etype).overlay(['l','m','t']).layout('Orb').cols(1)
Out[23]:
In [24]:
# Check also cross-secions
XS = np.abs(xrAF.XSrescaled)  # Should be in MB
# XS = np.abs(xrAF.XSraw)   # Unscaled/uncorrect XS
XS.name = 'XS/MB'
XS.hvplot.line(x='Eke').overlay('t').layout('Orb').cols(1)
Out[24]:

Versions

Code

In [25]:
# print(data.jobInfo['ePolyScat'][0])
# print('Run: ' + jobInfo['Starting'][0].split('at')[1])
In [26]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter', 'holoviews'])
Out[26]:
Fri Jan 07 13:47:48 2022 EST
OS Linux CPU(s) 64 Machine x86_64
Architecture 64bit Environment Jupyter
Python 3.7.10 (default, Feb 26 2021, 18:47:35) [GCC 7.3.0]
epsproc 1.3.1-dev xarray 0.15.0 jupyter Version unknown
holoviews 1.14.2 numpy 1.20.1 scipy 1.6.1
IPython 7.21.0 matplotlib 3.3.4 scooby 0.5.6
In [27]:
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1
* dev
  master
25c714d984e7f9deaf0bd39fd6171b79b185b592
In [28]:
# Check current remote commits
!git ls-remote --heads git://github.com/phockett/ePSproc
8fb518c63bfb4edc49ef994808f1c18fb68ff0c3	refs/heads/dependabot/pip/notes/envs/envs-versioned/pip-21.1
25c714d984e7f9deaf0bd39fd6171b79b185b592	refs/heads/dev
5278f0389878f6e7e5eb9b0c013ff3594c70938f	refs/heads/master
69cd89ce5bc0ad6d465a4bd8df6fba15d3fd1aee	refs/heads/numba-tests
ea30878c842f09d525fbf39fa269fa2302a13b57	refs/heads/revert-9-master

Notes

07/01/22 PH

  • Fixed/tested cos^2 values. (Maybe - now matching Reid 2018, but may not be correct in all cases?)
  • Updated plotting with HV for BLMs (in dev stage).
  • Checked XS values too - may have some -ve cases/phase issues to sort? (Cf. general AF code testing notes elsewhere.)

Versions

13/12/21 PH

STATUS: basics in place, need to trim code and fix a few things still.

Revisiting AF cases with (hopefully) better match to experimental alignment.

  • Add in cos^2(theta) metric.
  • Ballpark ADMs assuming Gaussian case.

TODO:

FIXED:

  • Note 'it' and 'labels' handling changed in recent versions of ePSproc, and degenerate dim hanling now correct!

09/06/21 PH

Quick look at initial ePS results for orb 11 (HOMO), 2.5eV step size.

TODO:

  • Currently AF code fails for Xarray = 0.17, issue with phaseCons setting to coord, sigh.
  • AF normalisation not correct here, due to doubly-degenerate state? Summing over 'it' not working correctly either, needs a debug (or missing some settings...?).
  • lmPlot() slice and cmapping to fix, in cases of missing/dropped dims, and for clims if possible.

For methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html

In [ ]: